FAIRE DES CARTOGRAMS DANS R

Dans cette séquance, nous proposons d’explorer les différentes méthodes disponibles dans R, pour réaliser des cartogrammes basés sur des données quantitatibes absolues.

Préalables

Installation des packages

install.packages(knitr)
install.packages(sf)
install.packages(mapsf)
install.packages(packcircles)
install.packages(cartogram)
install.packages(recmap)
install.packages(dplyr)
install.packages(cartogramR)
# install.packages("https://cran.r-project.org/src/contrib/Archive/cartogramR/cartogramR_1.0-1.tar.gz", repos = NULL, type = "source")

Environnement de travail

  1. Créez un nouveua projet
  2. Créez un script R
  3. Créez un répertoire data pour stocker vos données
  4. Créez un répertoire maps pour ranger vos cartes

Import et mise en forme des données

Les données utilisées ici sont des données INSSE sur la populatiion et les catégories socio professionelles en 2018 au niveau communal. L’espace d’étude est le département de l’Isère.

Le package sf (pour simple features), développé principalement par Edzer Pebesma (et Roger Bivan), permet de gérer les objets spatiaux dans R (projections, geotraoitepents, etc.). C’est le package qui a succedé au package sp.

library(sf)
communes <-
  st_read(
    "https://raw.githubusercontent.com/transcarto/rcartograms/main/data/isere.geojson",
    quiet = TRUE
  ) %>% st_transform(2154)
data <-
  read.csv(
    "https://raw.githubusercontent.com/transcarto/rcartograms/main/data/popisrere.csv",
    dec = ","
  )
communes = merge(x = communes[, c("id", "name", "geometry")],
                 y = data[, c("id",
                              "pop2018",
                              "agri",
                              "art",
                              "cadr",
                              "interm",
                              "emp",
                              "ouvr",
                              "retr")],
                 by = "id")
isere = st_union(communes)
knitr::kable(communes[c(0:10),], row.names = F, digits = 1)
id name pop2018 agri art cadr interm emp ouvr retr geometry
38001 Les Abrets en Dauphiné 6336 30.1 210.9 256.0 712.9 833.3 948.8 1397.2 MULTIPOLYGON (((900716.2 64…
38002 Les Adrets 1025 0.0 53.0 168.7 163.9 82.0 67.5 159.1 MULTIPOLYGON (((931354.5 64…
38003 Agnin 1127 0.0 28.9 72.2 168.6 139.7 134.9 192.6 MULTIPOLYGON (((844459.1 64…
38004 L’Albenc 1279 25.0 30.0 65.0 215.0 130.0 175.0 210.0 MULTIPOLYGON (((893395.8 64…
38005 Allemond 954 5.0 75.0 30.0 160.0 150.0 120.0 170.0 MULTIPOLYGON (((934142.3 64…
38006 Allevard 4062 0.0 176.7 323.2 535.2 469.6 469.6 984.2 MULTIPOLYGON (((948125.2 64…
38008 Ambel 28 10.0 0.0 0.0 0.0 0.0 0.0 10.0 MULTIPOLYGON (((930843.9 64…
38009 Anjou 1009 10.0 54.9 60.0 124.9 100.0 115.0 278.5 MULTIPOLYGON (((847739.9 64…
38010 Annoisin-Chatelans 686 14.8 44.4 54.2 113.3 64.1 78.9 113.3 MULTIPOLYGON (((875494.4 65…
38011 Anthon 1073 5.0 55.0 75.0 185.0 125.0 130.0 220.0 MULTIPOLYGON (((866538.7 65…

Exploration des représentations par packages : les classiques

Rappel : il s’agit ici de résoudre la question de la représentation de quantités brutes associées à des surfaces.

Rappel en carto “classique” avec le package mapsf

Le package mapsf permet de faire des cartes thématiques dans R. C’est le package qui succède au package cartography (qui n’est plus maintenu).

Trois solutions existent pour cartographier des quantités brutes associées à des surfaces en cartographie dite “classique” : la dot map, les points valués et les symboles proportionnels.

Création d’un template cartographique

library(mapsf)

col = "#c291bc"
credits = paste0("Bronner Anne-Christine & Nicolas Lambert, 2021\n",
                 "Source: IGN & INSEE, 2021")
theme = mf_theme(
  x = "default",
  bg = "#f0f0f0",
  tab = FALSE,
  pos = "center",
  line = 2,
  inner = FALSE,
  fg = col,
  mar = c(0, 0, 2, 0),
  cex = 1.9
)

template = function(title,
                    file,
                    note = "",
                    basemap = TRUE,
                    scale = TRUE) {
  mf_export(
    communes,
    export = "png",
    width = 1000,
    filename = file,
    res = 96,
    theme = theme,
    expandBB = c(-.02, 0, -.02, 0)
  )
  
  if (basemap == TRUE) {
    mf_shadow(
      x = communes,
      col = "grey50",
      cex = 1,
      add = TRUE
    )
    mf_map(
      communes,
      col = "#CCCCCC",
      border = "white",
      lwd = 0.5,
      add = TRUE
    )
  }
  
  mf_title(title)
  
  if (scale == TRUE) {
    mf_scale(
      size = 20,
      pos = "bottomright",
      lwd = 1.2,
      cex = 1,
      col = "#383838",
      unit = "km"
    )
  }
  
  mf_credits(
    txt = credits,
    pos = "bottomleft",
    col = "#1a2640",
    cex = 0.8,
    font = 3,
    bg = "#ffffff30"
  )
  
  if (note != "") {
    mf_annotation(
      x = c(885000, 6435000),
      txt = note,
      pos = "bottomleft",
      cex = 1.2,
      font = 2,
      halo = TRUE,
      s = 1.5
    )
  }
  
}
template("Template cartographique", "maps/template.png", note = "Département de\nl'Isère (38)")
mf_map(
  communes,
  col = col,
  border = "white",
  lwd = 0.5,
  add = TRUE
)
dev.off()

Dot density map ou carte par points

Il existe plusieurs solutions pour réaliser une carte par point dans R. Ici, nous vous proposeons de créer une fonction qui s’appuie sur st_sample() du package sf

dotdensitymap <- function(x,
                          var,
                          onedot = 1,
                          radius = 1) {
  x <- x[, c("id", var, "geometry")]
  x[, "v"] <- round(x[, var] %>% st_drop_geometry() / onedot, 0)
  dots <- st_sample(x, x$v, type = "random", exact = TRUE)
  circles <- st_buffer(dots, dist = radius)
  return (circles)
}
onedot = 500
dots = dotdensitymap(
  x = communes,
  var = "pop2018",
  onedot = onedot,
  radius = 300
)
template("Carte par points",
         "maps/dotdensity.png",
         note = paste0("Un point =\n", onedot, " habitants"))
mf_map(
  dots,
  col = col,
  border = "#520a2c",
  lwd = 0.5,
  add = TRUE
)
dev.off()

(Simili) Points Bertin ou points valués

Les points valués, appelés également “points Bertin” sont une variante de la dot map, avec deux-trois taille de points de référence pour s’adapter à de fortes disparités dans la distribution des effectifs. L’algorithme ci-après approche l’idée en distribuant les effectifs de chaque entité sur une grille sur des points de taille variable.

  grid = st_make_grid(communes, cellsize = 1000, square = TRUE)
  grid = st_sf(id = c(1:length(grid)), geometry = grid)
  ids = communes[, "id"] %>% st_drop_geometry()
  ctr = st_centroid(grid)
  
  for (i in 1:nrow(grid)) {
    x = st_within(
      x = ctr[i, ],
      y = communes,
      sparse = TRUE,
      prepared = TRUE
    )
    id = ids[unlist(x), ][1]
    if (length(id) == 0) {
      id = NA
    }
    grid[i, "id"] = id
  }
  
  grid <- grid[!is.na(grid$id), ]
  
  for (i in 1:nrow(grid)) {
    id = as.character(grid[i, "id"])[1]
    popAll = as.numeric(communes[communes$id == id, "pop2018"] %>% st_drop_geometry())
    count = nrow(grid[grid$id == id, ])
    grid[i, "pop"] = popAll / count
  }
  
  st_write(grid, "data/grid.geojson", delete_layer = TRUE)

La grille est maintenant disponible dans votre répertoire data.

Vous pouvez aussi aller chercher la grille précalulée ici

grid = st_read("https://raw.githubusercontent.com/transcarto/rcartograms/main/data/grid.geojson")
template("Points Bertin (grosso modo)", "maps/pointsbertin.png")
mf_map(communes, col = "#CCCCCC", border = "white", lwd = 0.5, add = FALSE)
mf_map(grid, var = "pop", col = col, border = "#520a2c", type = "prop",
       inches = 0.07, leg_title_cex = 1.2, leg_val_cex  = 0.8,
       leg_title = "Nombre d'habitants, 2018")

dev.off()

Symboles proportionnels

C’est la représentation la plus usuelle pour représenter des données quantitatives absolues.

template("Symboles proportionnels (mapsf)", "maps/prop.png")
mf_map(communes, var = "pop2018", col = col, border = "#6b4266", type = "prop",
       inches = 0.8, leg_title_cex = 1.2, leg_val_cex   = 0.8, symbol = "square",
       leg_title = "Nombre d'habitants, 2018")
dev.off()

template("Symboles proportionnels (mapsf)", "maps/prop2.png")
mf_map(communes, var = "pop2018", col = col, border = "#6b4266", type = "prop",
       inches = 0.8, leg_title_cex = 1.2, leg_val_cex   = 0.8,
       leg_title = "Nombre d'habitants, 2018")
dev.off()

Le cartogramme circulaire de Dorling avec le package packcircles

Le package packcirles propose 3 algorithmes simples pour déplacer des disques sur un plan 2D de telle sorte qu’ils ne se supperposent pas. Nous pouvons l’utiliser pour créer des cartogrammes de Dorling (Dorling 1996).

Création d’un ficher de données simplifié avec les coordonnées des centroïdes des communes.

dots = communes
st_geometry(dots) <- st_centroid(sf::st_geometry(dots),of_largest_polygon = TRUE)
dots <- data.frame(dots$id, dots["pop2018"], st_coordinates(dots))
dots = dots[,c("dots.id","X","Y","pop2018")]
colnames(dots) <- c("id","x","y","v")
dots <- dots[!is.na(dots$v),]
knitr::kable(dots[c(0:5),], row.names = F, digits = 1)
id x y v
38001 902026.9 6495943 6336
38002 934555.0 6466630 1025
38003 845708.9 6472468 1127
38004 892892.6 6460697 1279
38005 937029.6 6456880 954

La fonction circleRepelLayout() prend un ensemble de cercles dans un cadre de données et utilise la répulsion itérative pour essayer de trouver un arrangement sans chevauchement où tous les centres des cercles se trouvent à l’intérieur d’un rectangle de délimitation. Si aucun arrangement de ce type ne peut être trouvé dans le nombre maximum d’itérations spécifié, la dernière tentative est renvoyée.

library("packcircles")

k = 500 # pour ajuster la taille des cercles
itermax = 10 # nombre d'iterations

dat.init <- dots[,c("x","y","v")]
dat.init$v <- sqrt(dat.init$v * k)
simulation <- circleRepelLayout(x = dat.init, xysizecols = 1:3,
                                wrap = FALSE, sizetype = "radius",
                                maxiter = itermax, weights =1)$layout
knitr::kable(simulation[c(0:5),], row.names = F, digits = 1)
x y radius
902120.7 6496140 1779.9
934555.0 6466630 715.9
845708.9 6472468 750.7
892892.6 6460697 799.7
937029.6 6456880 690.7
circles <- st_buffer(sf::st_as_sf(simulation, coords =c('x', 'y'),
                      crs = sf::st_crs(communes)), dist = simulation$radius)
circles$v = dots$v

template("Dorling (packcircles)", "maps/dorling1.png", scale = FALSE)
mf_map(isere, col = "#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
mf_map(circles,col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()
## png 
##   2

Les 3 formes classiques avec le package cartogram

le package cartogram est développé par Sebastian Jeworutzki. Il propose trois méthodes : le cartogramme circulaire de Dorling, continu de Dougenik, Chrisman et Niemeyer et discontinu d’Olson).

library(cartogram)

Le cartogramme circulaire de Dorling

Proche d’un traitement graphique de l’information (le bubble chart), l’algorithme de Danny Dorling permet de visualiser les grandes masses de population en prenant en compte les positions relatives des entités les unes par rapport aux autres.

Dorling = cartogram_dorling(communes, "pop2018", k = 1.8)
template("Dorling (cartogram)", "maps/dorling2.png", scale = FALSE)
mf_map(isere, col = "#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
mf_map(Dorling, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

Le cartogramme discontinu d’Olson

Judith Olson, en réaction aux premiers cartogrammes continus de Tobler et notamment les problèmes de lisibilité, propose une variation proportionnelle de la surface tout en gardant sa forme “géographique” (Olson 1976)

Olson <- cartogram_ncont(communes, "pop2018", k = 1, inplace = TRUE)

inplace : Si VRAI, chaque polygone est redimensioné à sa place initiale, si FAUX, les multi-polygones sont centrés sur leur centroïde initial.

template("Olson (cartogram)", "maps/olson.png", basemap = FALSE, scale = FALSE)
mf_map(isere, col = "#CCCCCC30", border = "white", lwd = 0.5, add = TRUE)
mf_map(Olson, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

Le cartogramme continu Dougenik

James A. Dougenik, Nicholas R. Chrisman et Duane R.Niemeyer, sur la base des premiers travaux de W. Tobler, développent en 1985 un des deux algorithmes les plus implanté dans les outils (avec celui de Gastner et Newman développé en 2004). (Dougenik, Chrisman, and Niemeyer 1985)

x : un objet sf weight : om de la variable itermax : nombre d’itérations maximum si maxSizeError n’est pas atteint maxSizeError : la déformation s’arrete si l’erreur moyenne est inférieur à cette valeur prepare : mettre “adjust” permer d’acceler le temps de calcul. threshold : seuil pour la préparation des données

Dougenik <- cartogram_cont(communes, "pop2018", prepare = "none", itermax = 10, maxSizeError = 1.15)

Calcul des erreurs

sumarea = sum(as.numeric(st_area(Dougenik)))
sumpop = sum(Dougenik$pop2018)
Dougenik$error = (as.numeric(st_area(Dougenik)) /  sumarea) / (Dougenik$pop2018 / sumpop) * 100
summary(Dougenik$error)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   25.68   97.10  109.49  144.87  131.42 2731.97
bks = c(min(Dougenik$error),70,80,90,100,110,120,max(Dougenik$error))
cols = c("#d53e4f", "#f46d43","#fdae61","#fee08b","#e6f598","#abdda4", "#66c2a5")

Affichage de la carte

template("Dougenik (cartogram)", "maps/dougenik.png", basemap = FALSE, scale = FALSE)
mf_map(x = Dougenik, type = "choro",var = "error", pal = cols, breaks = bks, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

Les algos pour les cartogrammes continus avec le package cartogramR

library(cartogramR)

Le package cartogramR est développé par Pierre-Andre Cornillon et Florent Demoraes de l’université de Rennes. Le package propose les méthodes flow based cartogram (Gastner and Newman 2004) implanté dans Scapetoad, fast flow based cartogram (Gastner, Seguy, and More 2018) et rubber band based cartogram (Dougenik, Chrisman, and Niemeyer 1985). La méthode flow based cartogram est basée sur go_cart.

Attention, cartogramR n’est pas/plus sur le CRAN. L’archive est néanmoins disponible et permet l’installation.

la fonction precartogramR() aide à choisir la taille de la grille de déformation (plus elle est fine, plus le calcul sera long). On choisit donc a minima une grille telle que le minimum d’intersections est supérieur ou égal a un (ici un pas de grille de 256 peut faire l’affaire).

precarto <- precartogramR(communes, method = "GastnerSeguyMore")
summary(precarto)
##              L256      L512     L1024     L2048
## Min.      2.00000    5.0000   20.0000    82.000
## 1st Qu.  13.00000   52.0000  208.7500   837.000
## Median   19.00000   77.0000  306.0000  1226.500
## Mean     26.16992  104.6953  418.6816  1674.701
## 3rd Qu.  29.00000  115.2500  460.7500  1850.250
## Max.    405.00000 1613.0000 6469.0000 25884.000

la fonction cartogramR() permet de calculer le cartogramme en tant que tel selon les 3 méthodes proposées : GastnerSeguyMore" (ou “gsm),”GastnerNewman" (ou “gn),”DougenikChrismanNiemeyer" (ou “dcn”). L’option L=256 permet de choisir la taille de la grille. L’option grid=TRUE (pour les méthodes “gsm” et “gn”) permet d’afficher la grille de déformation. L’option maxit=50 permet de définir le nombre d’itérations max (défaut = 50).

Test de l’algorithme de Gastener et Newman de 2004, équivalent à celui implanté dans la solution java scapetoad.

GastnerSeguy <- cartogramR(communes, count="pop2018", method="GastnerSeguyMore", options=list(L=256, grid=TRUE, maxit = 5))
## WARNING criterion: 26.455251 > Objective: 0.010000
##  Increase maxit or decrease Objective

La fonction make_layer() permet de récupérer la grille de calcul (mais aussi la grille d’origine, les centroides, etc)

grid <- make_layer(GastnerSeguy, type = c("final_graticule"))

L’affichage des résultats peut se réaliser avec plot (via sf) ou, comme précédemment, avec le package mapsf.

template("Gastner, Seguy & More", "maps/gastnerseguy.png", basemap = FALSE, scale = FALSE)
mf_map(GastnerSeguy$cartogram, col = col, border = "white", lwd = 1, add = TRUE)
mf_map(grid, col = NA, border = "#6b4266", lwd = 0.05, add = TRUE)
dev.off()

La fonction residuals() permet de calculer les erreurs liées à la déformation (erreur relative : taille finale / taille theorique * 100)

table = cbind(GastnerSeguy$initial_data[,c("id","pop2018")] %>% st_drop_geometry(),
      orig_area = GastnerSeguy$orig_area,
      final_area = GastnerSeguy$final_area,
      errors = residuals(GastnerSeguy, type = "relative error")*100
      )
knitr::kable(table[c(0:10),], row.names = F, digits = 1)
id pop2018 orig_area final_area errors
38001 6336 27223532 39376361.8 -0.2
38002 1025 16353553 5636619.1 -11.7
38003 1127 8082831 7027719.4 0.2
38004 1279 10663444 7921991.8 -0.5
38005 954 56515131 5540617.6 -6.7
38006 4062 34980514 24517963.6 -3.0
38008 28 4789504 95852.1 -45.0
38009 1009 5121670 6220553.5 -1.0
38010 686 13149817 4281634.4 0.3
38011 1073 8596299 6774560.2 1.4

Export au format sf

st_write(as.sf(GastnerSeguy),"data/GastnerSeguy.geojson", delete_layer=TRUE)

Variations : encore plus de cartogrammes

Dots Cartograms

La méthode des dots cartograms est une représentation cartographique à l’intersection des cartogrammes de Dorling et de la carte par points. Les premières cartes utilisant cette méthodes ont été réalisées sur les données du Covid (voir et voir). Article à paraitre. Voir aussi sur Observable

dotcartogram = function(x,var,itermax,onedot,radius){
crs = sf::st_crs(x)
coords <- st_coordinates(st_centroid(sf::st_geometry(x),of_largest_polygon = TRUE))
x <- x[c("id",var)] %>% st_drop_geometry()
x <- data.frame(x, coords)
colnames(x) <- c("id","v","x","y")
x$v <- round(x$v/onedot,0)
x <- x[x$v > 0,]
dots <- x[x$v == 1,c("x","y","v")]
rest <-  x[x$v  > 1,c("x","y","v")]

nb <- nrow(rest)
  for (i in 1:nb){
    new <- rest[i,]
    new$v <- 1
    for (j in 1:rest$v[i]){ dots <- rbind(dots,new)}
  }
  
dots$x <- jitter(dots$x)
dots$y <- jitter(dots$y)
dots$v <- radius
  
simulation <- circleRepelLayout(x = dots, xysizecols = 1:3,
                                  wrap = FALSE, sizetype = "radius",
                                  maxiter = itermax, weights =1)$layout

circles <- st_buffer(sf::st_as_sf(simulation, coords =c('x', 'y'),
                                    crs = crs), dist = radius) 
return(circles)
}
onedot = 1000
dc = dotcartogram(x = communes, var = "pop2018", itermax = 120,
                  onedot = onedot, radius = 600)
template("Dots Cartogram", "maps/dotcartogram.png", note = paste0("Un point =\n",onedot," personnes"), scale = FALSE)
mf_map(isere, col = "#CCCCCC", border = "white", lwd = 0.5, add = TRUE)
mf_map(dc, col = col, border = "#6b4266", lwd = 0.8, add = TRUE)
dev.off()

Mosaic cartogram

Ici, on essaie de reproduire les cartogrammes de type mosaïque en trichant un peu… La forme choisie est l’hexagone. On est sur un hexagonal cartogram [voir]

On récupère un fond transformé par un algorithme de cartogramme continu

cartogram = st_read("https://raw.githubusercontent.com/transcarto/rcartograms/main/data/GastnerSeguy.geojson")
# Ou en local
# cartogram = st_read("data/GastnerSeguy.geojson")

Création d’une grille hexagonale, récupération des identificateurs des entités

grid = st_make_grid(cartogram, cellsize = 2000, square = FALSE)
grid = st_sf(id = rep("",length(grid)), geometry = grid)
ctr = st_centroid(grid)
ids = cartogram[,"id"] %>% st_drop_geometry()

for(i in 1:nrow(grid)){
  x = st_within(x = ctr[i,], y = cartogram, sparse = TRUE, prepared = TRUE)
  id = ids[unlist(x),][1]
  if (length(id) == 0){id = NA}
  grid[i,"id"] = id
}
grid <- grid[!is.na(grid$id),]

gridcartogram <- aggregate(x = grid, 
               by = list(grid$id),
               FUN = min)

Estimation de la valeur de chaque hexagone

varmax = sum(cartogram$pop2018)
nbcell = nrow(grid)
valcell = round(varmax / nbcell)

Création de la carte

template("Template cartographique", "maps/hexcartogram.png", note = paste0("Un hexagone ≈\n",valcell," habitants"),basemap = FALSE, scale = FALSE)
mf_shadow(x = grid, col = "grey50", cex = 1, add = TRUE)
mf_map(grid, col = col, border = "white", lwd = 0.5, add = TRUE)
mf_map(gridcartogram, col = NA, border = "#6b4266", lwd = 1, add = TRUE)
dev.off()

Hybride : tessalation et cartograme continu

NB : calculer les polygones de voronoi en R est un peu tricky. Voir et voir.

c <- st_centroid(communes)
v <- st_voronoi(x = st_union(c))
v <- st_intersection(st_cast(v), st_union(communes))
v <- st_join(x = st_sf(v), y = c, join=st_intersects)
v <- st_cast(v, "MULTIPOLYGON")
template("Polygones de voronoi", "maps/voronoi.png", basemap = FALSE, scale = FALSE)
mf_map(x = v, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

Ensuite, on peut calculer un cartogramme de Dougenik sur ces nouvelles géométries.

VDougenik <- cartogram_cont(v, "pop2018", prepare = "none", itermax = 15, maxSizeError = 1.15)
template("Dougenik & voronoi", "maps/dvoronoi.png", basemap = FALSE, scale = FALSE)
mf_map(x = VDougenik, col = col, border = "#6b4266", lwd = 0.5, add = TRUE)
dev.off()

L’impossibilité de réaliser un cartogramme rectangulaire avec le package recMap

A ce jour, il n’est pas possible d’automatiser la réalisation de cartogrammes de Demers ou rectangulaire, qui puissent construire et assembler un puzzle respectant surface, contiguité tout en recréant une forme globale à partir de formes carrées ou rectangulaires proportionnelles à des quantités brutes (stock).

Christian Panse a développé une première approche algorithmique que l’on trouve dans le package Recmap (Heilmann et al. 2004) (Panse 2016). L’algorithme va convertir les géométries des unités spatiales en rectangles dont la surface est définie par la variable thématique (stock). L’algorithme de RecMap est disponible ici. Développé en C++11, le package dépend des packages GA (>= 3.1), Rcpp (>= 1.0), sp(>= 1.3)

library(recmap)

Préparation des données

Création d’une fonction pour créer un fichier conforme à recmap

sfc_as_cols <- function(x, geometry, names = c("x","y")) {
  if (missing(geometry)) {
    geometry <- sf::st_geometry(x)
  } else {
    geometry <- rlang::eval_tidy(enquo(geometry), x)
  }
  stopifnot(inherits(x,"sf") && inherits(geometry,"sfc_POINT"))
  ret <- sf::st_coordinates(geometry)
  ret <- tibble::as_tibble(ret)
  stopifnot(length(names) == ncol(ret))
  x <- x[ , !names(x) %in% names]
  ret <- setNames(ret,names)
  dplyr::bind_cols(x,ret)
}

Extraction des centroïdes long-lat sous forme d’une table XY

coords <- st_centroid(communes)
coords <- st_transform(coords, crs="+proj=longlat +datum=WGS84 +ellps=WGS84")

Création du fichier pour recmap

df_recmap <- sfc_as_cols(coords)

Calcul des rectangles

rec_cartogram <- data.frame (x= df_recmap$x,
                             y=df_recmap$y,
                             # make the rectangles overlapping by correcting lines of longitude distance
                             dx = sqrt(df_recmap$pop2018) / 90 / abs((0.65 * 60 * cos(df_recmap$y*pi/180))),
                             dy = sqrt(df_recmap$pop2018) / 90 / (0.65 * 60),
                             z = sqrt(df_recmap$pop2018),
                             name = df_recmap$id)

Création de la carte

template("RecMap", "maps/recmap1.png", basemap = FALSE, scale = FALSE)
plot.recmap(rec_cartogram, col = NA, border = col, lwd=4,  col.text = col)
dev.off()

# cartog <- recmap(df)
# template("RecMap", "maps/recmap2.png", basemap = FALSE, scale = FALSE)
# plot(cartog[!cartog$name %in% c("IS","MT"),],  col = col, border = "white")
# dev.off()

Attention, il ne fonctionne que sur un nombre limité d’unités territoriales. Ici, l’algo ets incapable de déplacer les rectangles pour creer un cartogramme. Seule solution : exporter en svg et les déplacer à la main.

Avec un nombre plus réduit d’unités, il est possible de mener le process jusqu’au bout. Essayons sur l’Europe…

europe <-
  st_read(
    "https://raw.githubusercontent.com/transcarto/rcartograms/main/data/europe.geojson"
  ) %>% st_transform(3035)
coords = data.frame(st_coordinates(st_centroid(st_geometry(europe))))
bb <- lapply(st_geometry(europe), function(x) {
  st_bbox(x)
})
dx <- unlist(lapply(bb, function(x) {
  x[3] - x[1]
})) / 2
dy <- unlist(lapply(bb, function(x) {
  x[4] - x[2]
})) / 2

df <- data.frame(
  x = coords$X,
  y = coords$Y,
  dx = dx,
  dy = dy,
  z = europe$pop2008,
  name = europe$id
)
template("", "maps/recmap2.png", basemap = FALSE, scale = FALSE)
plot.recmap(
  df,
  col = NA,
  border = col,
  lwd = 4,
  col.text = col
)
dev.off()

cartog <- recmap(df)
template("", "maps/recmap3.png", basemap = FALSE, scale = FALSE)
plot(cartog[!cartog$name %in% c("IS", "MT"),],  col = col, border = "white")
dev.off()

NB : IS et MT n’ont pu être placés. On les supprime à l’affichage.

Le cartogramme comme système de projection

S’il résoud la question de la variation de taille en implantation surfacique, le cartogramme permet également de “corriger” le fond de carte chroplèthe lorsque l’on cartographe des phénomènes socio-démographiques (cartographie électorale, données insee, de santé…).

Il suffit d’appeler lors de la cartographie de la variable attributaire le fond transformé. A priori, le fond transformé correspond à la valeur du dénominateur du taux cartographié (le nombre d’inscrits sur les listes électorales pour cartographier le taux d’abstention, la population active pour un taux de chômage, etc.). Néanmoins lors d’étude sur plusieurs variables socio-démo, on se base souvent sur un fond unique de référence (en général la population résidentielle).

Exemple cartographie de la part de cadres et d’ouvriers

Du coup là on verra pour les variables qu’on propose, juste faire gaffe d’avoir le bon dénominateur pour le calcul du taux, cf. population 18-65 pour calculer un taux à partir de csp.

Comparaison taux de cadres et ouvriers

Calcul des taux et discrétisation

communes$ouvr_tx <- communes$ouvr / communes$pop2018 * 100
communes$cadr_tx <- communes$cadr / communes$pop2018 * 100
bks <- c(0, 0.5, 5, 10, 15, 20, 25, 100)
cols <-
  c("#ffffff",
    "#e5f5f9",
    "#ccece6",
    "#66c2a4",
    "#238b45",
    "#006d2c",
    "#00441b")

71% de cadres à Oulles, 6 habitants ;-)

Calcul des fonds de carte

Dougenik <- cartogram_cont(communes, "pop2018", prepare = "none", itermax = 10, maxSizeError = 1.15)
Olson <- cartogram_ncont(communes, "pop2018", k = 1, inplace = TRUE)
Dorling = cartogram_dorling(communes, "pop2018", k = 1.8)

Cartographies

Carte choroplèthe = quel fond de carte choisir pour représenter un taux ? géographique, cartogramme continu, discontinu ou de Dorling

Pour les ouvriers

par(mfrow = c(2, 2), mar = c(0,0,0,0))
mf_map(x=communes,var = "ouvr_tx", type = "choro", pal = cols, breaks = bks)
mf_map(x=Dougenik,var = "ouvr_tx", type = "choro", pal=cols, breaks = bks)
mf_map(x=Dorling,var = "ouvr_tx", type = "choro",  pal=cols, breaks = bks)
mf_map(x=communes, col=NA)
mf_map(x=Olson,var = "ouvr_tx", type = "choro", pal=cols, breaks = bks, add = TRUE)

Pour les cadres

par(mfrow = c(2, 2), mar = c(0,0,0,0))
mf_map(x=communes,var = "cadr_tx", type = "choro", pal = cols, breaks = bks)
mf_map(x=Dougenik,var = "cadr_tx", type = "choro", pal=cols, breaks = bks)
mf_map(x=Dorling,var = "cadr_tx", type = "choro",  pal=cols, breaks = bks)
mf_map(x=communes, col=NA)
mf_map(x=Olson,var = "cadr_tx", type = "choro", pal=cols, breaks = bks, add = TRUE)

Comparaison cadres-ouvriers sur fond transformé (peut-être faut faire versus fond géo)

par(mfrow = c(1, 2), mar = c(0,0,0,0))
mf_map(x=Dougenik,var = "cadr_tx", type = "choro", pal=cols, breaks = bks)
mf_map(x=Dougenik,var = "ouvr_tx", type = "choro", pal=cols, breaks = bks)

Versus comparaison par la forme qui reproduit en général la répartition de la population lorsque l’on est sur des données démographiques. = cartogramme directement sur les effectifs.

Changement d’échelle !

Cartogramme et échelle spatiale sont fondamentalement liés. La lecture de toute carte, comme du cartogramme est plus ou moins liée à la familiarité avec l’espace représenté. Dans le cas du cartogramme, l’étendue de l’espace cartographié, le niveaux d’agrégation et la méthode choisie permettent de chercher une solution adaptée au phénomène à représenter. Le cartogramme rectangulaire ou de Demers est probablement plus adaptés lorsque le nombre d’entité à cartographier est limité, alors que les cercles de Dorling s’adaptent à toute situation.

Cartogramme et étendue spatiale

Explorez les cartogrammes à l’échelle mondiale

Quel fond de carte pour la représentation du monde ?

world <- st_read("https://raw.githubusercontent.com/transcarto/rcartograms/main/data/world_countries.geojson", quiet = TRUE ) %>%
  st_transform( "+proj=bertin1953")

world <- world[world$ISO3 != "ATA",]
knitr::kable(world[c(0:10),], row.names = F, digits = 1)
ISO2 ISO3 ISONUM NAMEen NAMEfr UNRegion GrowthRate PopDensity PopTotal JamesBond geometry
GQ GNQ 226 Equatorial Guinea Guinée-Équatoriale Africa 2.9 30.1 845060 0 MULTIPOLYGON (((-514804.2 -…
TV TUV 798 Tuvalu Tuvalu Oceania 0.3 330.5 9916 0 MULTIPOLYGON (((11846963 52…
MV MDV 462 Maldives Maldives Asia 1.7 1212.2 363657 0 MULTIPOLYGON (((5515581 -20…
AD AND 20 Andorra Andorre Europe -1.9 149.9 70473 0 MULTIPOLYGON (((-1013908 15…
AG ATG 28 Antigua and Barbuda Antigua-et-Barbuda America 1.0 208.7 91818 0 MULTIPOLYGON (((-6420562 59…
AT AUT 40 Austria Autriche Europe 0.3 103.7 8544586 4 MULTIPOLYGON (((27563.98 73…
BE BEL 56 Belgium Belgique Europe 0.6 373.2 11299192 0 MULTIPOLYGON (((-621349.5 1…
BH BHR 48 Bahrain Bahren Asia 1.4 1812.2 1377237 0 MULTIPOLYGON (((2804520 -10…
BS BHS 44 Bahamas Bahamas America 1.2 38.8 388019 3 MULTIPOLYGON (((-6563129 25…
BB BRB 52 Barbados Barbade America 0.3 661.0 284215 0 MULTIPOLYGON (((-6530926 70…

Un premier point de vue sur le monde (projection Bertin, par Fil)

col = "#c291bc"
credits = "Vous, 2021"
theme = mf_theme(x = "default", bg = "#f0f0f0", tab = FALSE, 
                   pos = "center", line = 2, inner = FALSE, 
                   fg = col, mar = c(0,0, 2, 0),cex = 1.9)

templateworld = function(title, file){

  mf_export(
    world,
    export = "png",
    width = 1000,
    filename = file,
    res = 96,
    theme = theme, 
    expandBB = c(-.02,0,-.02,0)
)

  mf_title(title)

  mf_credits(
    txt = credits,
    pos = "bottomleft",
    col = "#1a2640",
    cex = 0.8,
    font = 3,
    bg = "#ffffff30"
  )
}

Votre carte de base

templateworld("Le monde en projection Bertin (thx Fil)", "maps/world.png")
mf_map(world, col =col, border = "white", lwd = 0.5, add = TRUE)
dev.off()

Lorsque l’on travaille avec les cartogrammes, les messages d’erreurs les plus courants sont liés à des valeurs absentes, des polygones mal fermés ou des contiguités imparfaites. Cela se discute, mais on peut considérer qu’il est plus “propre” de supprimer les entités avec des valeurs absentes ou nulles du jeu de données : en effet, elles ne sont pas concernées par le phénomène représenté.

world_sansNA = na.omit(world)

Dougenick

templateworld("La population mondiale", "maps/worldDougenick.png")
mf_map(x=cartogram_cont(world_sansNA, "PopTotal"), col =col, border = "white", lwd = 0.5, add = TRUE)
dev.off()

Dorling

templateworld("La population mondiale", "maps/worldDorling.png")
mf_map(cartogram_dorling(world_sansNA, "PopTotal"), col =col, border = "white", lwd = 0.5, add = TRUE)
dev.off()

Olson

templateworld("La population mondiale", "maps/worldOlson.png")
mf_map(x=cartogram_ncont(world_sansNA, "PopTotal"), col =col, border = "white", lwd = 0.5, add = TRUE)
dev.off()

Explorez les cartogrammes à l’échelle européenne

europe <- st_read("https://raw.githubusercontent.com/transcarto/rcartograms/main/data/europe.geojson") %>% st_transform(3035)

Les formes du monde (le fond de carte en question)

par(mfrow = c(2, 2), mar = c(0,0,0,0))
mf_map(europe, col = col)
mf_map(x=cartogram_cont(europe, "pop2008"), col = col)
mf_map(x=cartogram_dorling(europe, "pop2008"), col = col)
mf_map(x=cartogram_ncont(europe, "pop2008"), col = col)

dev.off()

A vous de jouer…

Et pour pour les afficionados de R, quelques pistes à creuser : * travailler l’habillage de la carte : comment automatiser l’ajout d’une légende (surface de référence) * améliorer la qualité du cartogramme continu (cf. mesure d’erreur) : repasser l’algorithme sur le fond transformé se heurte souvent à des erreurs de topologie créées par la transformation. Ceci nécessite de passer l’équivalent du “réparer les topologies” de QGIS sur les fonds transformés (opération a éventuellement effectuer par défaut pour un fond plus “propre”).

Un mot sur l’anticartogramme

Eduardo Dutenkefer a exploré l’esthétique des cartogrammes et notamment imaginé l’anti-cartogramme, qui inverse le poids des entités : la surface varie en proportion inverse à la valeur cartographiée.

Rappel. Le cartogramme a pour objectif par la variation de taille, qui a un sens bien défini en cartographie, basé sur la perception visuelle (la surface permet d’estimer la part de chaque ou d’un groupe d’entités sur le total), de révéler les inégalités de répartition sur un territoire.

Si on s’en tient au traitement cartogramme de la donnée, l’anticartogramme ne sert à rien dans le cadre d’une analyse des inégalités territoriales. Il serait faux de dire que l’anticartogramme représenterait le contraire (excepté dans le cas où deux variables représenterait la part d’un tout : p. ex. les votes exprimés entre deux candidats, chômeurs et actifs dans la population active, situation où nous sommes de facto en présence d’une relation inverse entre les deux variables).

Représenter le non poids des lieux ? Un cartogramme de la non population ne représente pas nécessairement les zones dites rurales. A l’échelle mondiale, l’anti-cartogramme des riches n’est pas le cartogramme des pauvres. Les inégalités de répartition de la non-population, du non-PIB, du non CO2, de la même manière que la visualisation des non flux ou une carte des anti flux, où le flux le plus insignifiant devient le plus massif. Pour révéler et analyser les zones où les effectifs sont faibles, les stocks peu élevés, le traitement de la donnée s’effectue par le filtrage.

Visualiser le “non poids,” créer la “non image” d’un espace peut soutenir un discours artistique, métaphorique sur un espace et non pas à cartographier une variable quantitative.

(Poncet 2011) Je mettrais plutôt la thèse de Dutenkefer ou son article avec Padovesi Fonseca antérieures au topo de Poncet. Si j’avais eu un peu le temps je les aurais contactés car ils auraient sûrement été très intéressés par Transcarto.

Calcul de la valeur du “non phénomène” cartographié (1/)

communes$inverse = 1/communes$pop2018

Bibliographie

Dorling, D. 1996. “Area Cartograms: Their Use and Creation, Concepts and Techniques in Modern Geography.”
Dougenik, James A, Nicholas R Chrisman, and Duane R Niemeyer. 1985. “An Algorithm to Construct Continuous Area Cartograms.” The Professional Geographer 37 (1): 75–81.
Gastner, Michael T, and Mark EJ Newman. 2004. “Diffusion-Based Method for Producing Density-Equalizing Maps.” Proceedings of the National Academy of Sciences 101 (20): 7499–7504.
Gastner, Michael T, Vivien Seguy, and Pratyush More. 2018. “Fast Flow-Based Algorithm for Creating Density-Equalizing Map Projections.” Proceedings of the National Academy of Sciences 115 (10): E2156–64.
Heilmann, Roland, Daniel A Keim, Christian Panse, and Mike Sips. 2004. “Recmap: Rectangular Map Approximations.” In IEEE Symposium on Information Visualization, 33–40. IEEE.
Olson, Judy M. 1976. “Noncontiguous Area Cartograms.” The Professional Geographer 28 (4): 371–80.
Panse, Christian. 2016. “Rectangular Statistical Cartograms in r: The Recmap Package.” arXiv Preprint arXiv:1606.00464.
Poncet, Patrick. 2011. “Antigéographie. L’anticartogramme Et Ses Interprétations.”